function [x_k, stats] = L_BFGS(x_k, fun, options, plot_fun)
% [x_k, stats] = L_BFGS(x_k, fun, options, plot_fun)
% 
% L_BFGS is an implementation of the large scale BFGS algorithm [1, p.
% 226]. Line search is performed as Armijo Backtracking [1, p. 41]. 
%
% [1] J. Nocedal and S. J. Wright, 
%     Numerical Optimization, 
%     New York: Springer, 2000
%
% See also: Backtracking, Check_convergence
%
% 2013-07-25 --------------------------------------------------------------
% 
% Matthias J. Ehrhardt
% CMIC, University College London, UK
% 
% matthias.ehrhardt.11@ucl.ac.uk
% 
% -------------------------------------------------------------------------

    % treat the options
    if ~exist('options', 'var'); options = struct; end;
        
    if ~isfield(options, 'x_tolerance'),    options.x_tolerance     = 1e-6;     end; % Choose as -Inf to switch this off
    if ~isfield(options, 'n_iter'),         options.n_iter          = 20;       end;
    if ~isfield(options, 'n_dt_iter'),      options.n_dt_iter       = 20;       end;
    if ~isfield(options, 'm'),              options.m               = 3;       end;
    if ~isfield(options, 'delta'),          options.delta           = 1;        end;
    if ~isfield(options, 'verbose'),        options.verbose         = false;    end;
    
    n_iter  = options.n_iter;
    m       = options.m;
    verbose = options.verbose;
    
    [fun_k, grad_k, norm_grad_k, extra_k] = fun(x_k);
    max_dx = 0;
    
    % for statistics
    if nargout == 2
        stats.fun               = zeros(n_iter+1, 1);
        stats.norm_grad         = zeros(n_iter+1, 1);
        stats.extra             = cell(n_iter+1, 1);          
        stats.dt                = zeros(n_iter, 1);
        stats.dx                = zeros(n_iter, 1);
        stats.time              = zeros(n_iter, 1);
        
        stats.converged         = -1;
        
        stats.fun(1)            = fun_k;
        stats.norm_grad(1)      = norm_grad_k;
        stats.extra{1}          = extra_k;       
    end
    
    if nargout == 2 || options.verbose
        timer = tic;
    end
    
    % auxiliary variables for L_BFGS
    s       = cell(n_iter, 1);
    y       = cell(n_iter, 1);
    rho     = cell(n_iter, 1);
        
    % display progress
    if verbose
        display(' --- start L-BFGS iterations');
        fprintf('%s\n', Output_string(0, toc(timer), fun_k, 0, norm_grad_k, 0)); 
    end
    
    if nargin == 4
        plot_fun(0, x_k)
    end
    
    for k = 1 : n_iter        
        % Choose H_k0 (for example by using equation 9.6)
        if k > 1
            gamma_k = (s{k-1}'*y{k-1})/(y{k-1}'*y{k-1});
        else
            gamma_k = options.delta / norm_grad_k;
        end
                
        % Compute p_k via algorithm 9.1, p.225
        p_k = - L_BFGS_two_loop_recursion(k, m, gamma_k, grad_k, s, y, rho);
        
        [dt_k, x_kk, fun_kk, grad_kk, norm_grad_kk, extra_kk] = ...
            Backtracking(options.n_dt_iter, x_k, fun, fun_k, ...
            p_k, false);

        % Converged? (because of a small stepsize)
        if isnan(dt_k)
            verbose = true;
            if verbose; fprintf('   >>> Converged! (because no stepsize found which would reduce the function value)\n'); end;
            if nargout == 2; stats.converged = k-1; end;
            return
        end
        
        % Converged?
        converged = Check_convergence(k, x_k, x_kk, options, true);
        
        if converged
            if nargout == 2; stats.converged = k-1; end;
            return;
        end
        
        % If k is large enough, clear unused space
        if k > m
            s{k-m} = [];
            y{k-m} = [];
        end
       
        % Compute and safe new values for s_k and y_k
        s{k} = x_kk - x_k;
        y{k} = grad_kk - grad_k;
        % and for computational reasons rho_k (that differs from the book)
        rho{k} = 1./ (y{k}' * s{k});
        
        if nargin == 4
            plot_fun(k, x_kk);
        end
              
        % update for next iteration
        max_dx = norm(x_kk - x_k, 'inf');
        x_k = x_kk;
        grad_k = grad_kk;
        norm_grad_k = norm_grad_kk;
        fun_k = fun_kk;
        
        % display progress
        if verbose; fprintf('%s\n', Output_string(k, toc(timer), fun_kk, dt_k, norm_grad_kk, max_dx)); end;

        % stats
        if nargout == 2
            stats.fun(k+1)           = fun_kk;
            stats.norm_grad(k+1)     = norm_grad_kk;
            stats.extra{k+1}         = extra_kk;
            stats.dt(k)              = dt_k;
            stats.dx(k)              = max_dx;
            stats.time(k)            = toc(timer);
        end
    end
    
    % display progress
    verbose = true; 
    if verbose; 
        fprintf('%s\n', Output_string(k, toc(timer), fun_kk, dt_k, norm_grad_kk, max_dx));
        warning('L_BFGS: NOT converged because the maximum number of iterations reached.\n');
    end;
end

% result is H_k * grad_x_k
function r = L_BFGS_two_loop_recursion(k, m, gamma_k, grad_x_k, s, y, rho)
    alpha = zeros(m,1);
    q = grad_x_k;
    for i = k-1 : -1 : k-m
        if i > 0            
            alpha(i) = rho{i} * s{i}'*q;
            q = q - alpha(i)*y{i};
        end
    end
    
    r = gamma_k * q;
    
    for i = k-m : k-1
        if i > 0
            beta = rho{i} * y{i}'*r;
            r = r + s{i} * (alpha(i) - beta);
        end
    end
end